双重差分法（DID）安慰剂检验的做法：随机抽取500次?

“安慰剂”（placebo）一词来自医学上的随机实验，比如要检验某种新药的疗效。此时，可将参加实验的人群随机分为两组，其中一组为实验组，服用真药；而另一组为控制组，服用安慰剂（比如，无用的糖丸），并且不让参与者知道自己服用的究竟是真药还是安慰剂，以避免由于主观心理作用而影响实验效果。

双重差分法（DID）安慰剂检验的核心思想就是虚构处理组或者虚构政策时间进行估计，如果虚构情况下“伪政策虚拟变量”的系数依然显著，那么就说明原来的估计结果很有可能出现了偏误，我们的被解释变量y的变动很有可能是受到了其他政策或者随机性因素的影响。

说到虚构，那么自然是可以随机虚构，也可以不随机虚构（作者自己设定）。当然，我更推荐的还是随机虚构处理组或者是政策时间的方法。由于我们使用的数据基本都是“大N小T”型的短面板数据，所以随机虚构政策时间没什么意义，文献一般做法都是将政策年份统一提前2年或3年重新进行回归，看看政策虚拟变量系数是否依然显著。我们更多地还是随机虚构处理组，具体做法就是随机选取个体作为处理组，重复500次或者1000次，看看“伪政策虚拟变量”的系数是否显著。

数据来源

石大千等（2018）发表在《中国工业经济》的论文《智慧城市建设能否降低环境污染》使用DID方法评估了智慧城市建设对城市环境污染的影响，《中国工业经济》期刊官网公布了这篇论文使用的数据和代码。接下来，我就使用这篇论文的数据，给大家分享一下双重差分法（DID）安慰剂检验中随机虚构处理组这种方法的Stata操作。

原文信息

石大千,丁海,卫平,刘建江.智慧城市建设能否降低环境污染[J].中国工业经济,2018(06):117-135.

随机虚构处理组的Stata操作

双重差分法（DID）安慰剂检验的一般做法就是随机选取个体作为处理组，重复500次或者1000次，看看“伪政策虚拟变量”的系数是否显著。在石大千等（2018）这篇论文中，处理组有32个城市，控制组有165个城市，所以我们需要从197个城市中随机选取32个城市作为“伪处理组”，假设这32个城市是智慧城市试点，其他城市为控制组，然后生成“伪政策虚拟变量”（交互项）进行回归。

重复进行500次，我们就会得到500次回归结果（包含“伪政策虚拟变量”估计系数、标准误和p值），最后我们可以绘制出500个“伪政策虚拟变量”估计系数的分布及相应的p值，直观展示安慰剂检验的结果。

说起来很简单，但操作起来还是蛮困难的，我们要解决如下两个关键问题：

（1）如何从197个城市中随机选取32个城市作为“伪处理组”，然后生成“伪政策虚拟变量”？因为是面板数据，所以我们首先只保留一期数据，然后使用sample命令从中随机抽取32个样本，保留所抽取样本的id编号，与原数据进行匹配，最终匹配上的就是我们的“伪处理组”样本，未匹配上的就是控制组样本。

（2）如何存储500次回归中的估计系数、标准误和p值？在Stata中，我们可以生成三个500行1列的矩阵（矩阵中的元素初始值为0），来分别存储500个“伪政策虚拟变量”的估计系数、标准误和p值。每回归一次，就将估计系数、标准误和p值分别存储到三个矩阵的对应位置。

最终代码

*-安慰剂检验-虚构处理组
mat b = J(500,1,0) //* 系数矩阵
mat se = J(500,1,0) //* 标准误矩阵 
mat p = J(500,1,0) //* P值矩阵

*循环500次
forvalues i=1/500{
 use smart_city2018.dta, clear
 xtset id year  //面板数据声明
 keep if year==2005   //保留一期数据
 sample 32, count   //随机抽取32个城市
 keep id  //得到所抽取样本的id编号
 save match_id.dta, replace   //另存id编号数据
 merge 1:m id using smart_city2018.dta //与原数据匹配
 gen treat = (_merge == 3) //将所抽取样本赋值为1，其余为0，得到政策分组虚拟变量
 gen period = (year >= 2012) //生成政策时间虚拟变量
 gen did = treat*period
 reghdfe lnrso did $xlist ,absorb(id year) vce(cluster id) 

 * 将回归结果赋值到对应矩阵的对应位置
 mat b[`i',1] = _b[did]
 mat se[`i',1] = _se[did]

 * 计算P值并赋值于矩阵
 mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did]/_se[did]))
}

* 矩阵转化为向量
svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)
* 删除空值并添加标签
drop if pvalue1 == .
label var pvalue1 p值
label var coef1 估计系数
keep coef1 se1 pvalue1 

绘图代码

twoway (kdensity coef1) (scatter pvalue1 coef1, msymbol(smcircle_hollow) mcolor(blue)), ///
title("Placebo Test") ///
xlabel(-0.4(0.1)0.4) ylabel(,angle(0)) ///
xline(-0.171, lwidth(vthin) lp(shortdash)) xtitle("Coefficients") ///
yline(0.1,lwidth(vthin) lp(dash)) ytitle(p value) ///
legend(label(1 "kdensity of estimates") label( 2 "p value")) ///
plotregion(style(none)) ///无边框
graphregion(color(white)) //白底

上图展示了500个“伪政策虚拟变量”估计系数的分布及相应的p值，其中X轴表示“伪政策虚拟变量”估计系数的大小，Y轴表示密度值和p值大小，曲线是估计系数的核密度分布，蓝色圆点是估计系数对应的p值，垂直虚线是DID模型真实估计值-0.171，水平虚线是显著性水平0.1。

从图中可以看出，估计系数大都集中在零点附近，大多数估计值的p值都大于0.1（在10%的水平上不显著），这表明我们的估计结果不太可能是偶然得到的，因而不太可能受到了其他政策或者随机性因素的影响。不过，这个安慰剂检验的结果并不是那么理想，可能还需要在之前的一些环节做些微调，最好能让我们从图中看到DID模型真实估计值是明显的异常值。

这里，我只是提供了一种主流的安慰剂检验的做法，希望大家不要用“定式思维”来看待安慰剂检验。除此之外，还有其他各式各样的安慰剂检验做法，要根据政策情况、样本情况等方面进行具体分析。需要本篇推文使用的数据和代码的朋友，请在公众号后台对话框内回复关键词“placebo”。